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Abstract 

The problem of calculating multicanonical parameters recursively is discussed. I describe in 
detail a computational implementation which has worked reasonably well in practice. 
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1 Introduction 



Recently Monte Carlo (MC) sampling with respect to unconventional ensembles has received 
some attention [1-20] . In the multicanonical ensemble |j| one samples configurations such 
that exact reconstruction of canonical expectation values becomes feasible for a desired tem- 
perature range. This requires a broad energy distribution, and leaves innovative freedom 
concerning the optimal shape [[/J. Considerable practical experience exists only for the uni- 
form energy distribution, where one samples such that: 

(a) The energy density is flat in a desired range 

P(E) = const for E min < E < E max . (1) 

(b) Each configuration of fixed energy E appears with the same likelihood. 

It should be noted that condition (b) is non-trivial. A simple algorithm [2T| exists to 



achieve (a), but which give up (b). Exact connection to the canonical ensemble is then lost. 
Such algorithms are interesting for hard optimization problems, but unsuitable for canonical 
statistical physics. The present paper focuses on achieving (a) and (b). 

The average computer time r, measured in updates, which it takes to proceed from 
E min to E max and back has been named "tunneling time" 0, 0- It should be noted that 
the method overcomes free energy barriers actually not by a tunneling process, but through 
moving along valleys, which are connected to the disordered phase. Once an updating scheme 
is given, like standard Metropolis, it is an interesting theoretical question to find the weight 
factors which minimize the tunneling time. It is by no means clear that this will be the 
uniform choice (1), on which the present paper is focused. 

Multicanonical and related sampling has allowed considerable gains in situations with 
"supercritical" slowing down. Such are: 



(a) First order transitions (T|, ||, for a recent review see J20 . 
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(b) Systems with conflicting constraints, such as spin glasses |2|, ^, [T7|, [H| or proteins 

pi ph. 



To achieve a flat energy distribution, the appropriate unnormalized weight factor w(E) 
is the inverse spectral density w(E) = n'~ 1 {E) 1 just like the weight factor for canonical MC 
simulations is the Boltzmann factor w B (E) = exp(— f3E). Now, the spectral density is a- 
priori unknown. Otherwise we would have solved the problem in the first place. Presumably, 
reluctance about simulations with an a-priori unknown weight factor is the main reason why 



the earlier umbrella sampling |22fl never became popular in statistical physics. 



For first order phase transitions the problem of the a-priori unknown weight factor is 



rather elegantly overcome by means of finite size scaling (FSS) methods [|T], |9|, [TO], |T^, |20 
A sufficiently accurate estimate is obtained by extrapolation from the already simulated 
smaller lattices. The smallest lattices allow still for efficient canonical simulations. 

For systems with conflicting constraints the situation is less satisfactory. For instance 
for spin glasses one has to perform the additional average over quenched random variables 
(which are the exchange coupling constants). Different choices of these random variables 
define different realizations of the same system. For the Edward-Anderson Ising spin glass 
it turned out [0, |17]] that, even for identical lattice sizes, different realizations need different 
weight factors. Each system requires a new estimate of the weight factors with no a-priori 
information available. To achieve this, a recursion 23 was introduced by Celik and the author 
[0]. However, details of the recursion (see section 3) may need considerable attention by 
hand. This attention is possible when only a few lattices are simulated, but impractical 
when hundreds or even thousands of different realizations have to be handled. This renders 
it inconvenient for more complicated situations, like the 3d Edwards-Anderson Ising (EAI) 
spin glass. 

Consequently, the recursion actually used in Ref.[|17|, where multicanonical simulations 



were performed for more then 1,500 different realizations of the EAI model, differed from 
the one described in [E|. The main purpose of this article is to describe this particular 



approach. In each recursion step the statistical information from all previous runs is used 
directly for estimating the multicanonical parameters as well as for noise reduction. Further, 
the recursion turned out to be robust. Little attention by "hand" was needed. However, 
no claim is made that it is in any sense optimal (actually the author is considering various 
improvements). It is supposed to be a reasonable starting point to provide a running code 
quickly. 

The paper is organized as follows: In section 2 generalized Ising models and related pre- 
liminaries are introduced. Mainly for pedagogical reasons I focus on them for examples of 
this paper. It is clear that generalization to other systems is straightforward, although possi- 
bly tedious for continuous systems. In section 3 I introduce the multicanonical method, and 
discuss the recursions given in the literature 0, I3]| . Section 4 describes the recursion which 
I invented for the simulations of |]l~Tfl , and section 5 illustrates its performance. Summary 
and conclusions follow. The appendix gives and explains a corresponding program listing. 



2 Generalized Ising Models 

Let us consider a <i-dimensional hypercubic lattice of volume V = N = L d with periodic 
boundary conditions. Spins Sj = ±1 are located at the N sites, and exchange interactions 
Jik = ±1 at the dN links of the lattice. The energy of generalized Ising models is given by 

E = -J2 JikSiSj, (2) 

ik 

where the sum is over the nearest neighbors. For J^. = 1 the standard Ising ferromagnet (IF) 
is recovered. When the are quenched random variables, one obtains the EAI spin glass. 
I confine the subsequent discussion to these two situations, although there are other cases of 
interest [Q. Let us further restrict the EAI spin glass to the situation J2<ik> Jik = 0. 
The partition function may be written as 

Z((3) = J2n(E)e^ E } (3) 

E 
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where n(E) is the spectral density |[25| , more precisely the number of configurations (or 
states) with energy E. As the system has 2 N different states, this implies the normalization 

J2n(E) = 2 N . (A) 

E 

The lowest possible energy is —dN, reached when each link contributes JikSiSk = 1. For 
the IF this is achieved with either all spins up (+1) or all spins down (—1). For a generic 
configuration the possible energy increments under the flip of a single spin are 

AE = 0, ±4, ±Ad. (5) 

Consequently n(E) may take non-zero values for 

E = -dN, -dN + 4, 0, dN — 4, dN. (6) 

For instance for the IF n(-dN) = 2, n(-dN + 4) = 0, and n(-dN + Ad) = N. For a 
typical EAI spin glass configuration the groundstate energy E min is considerably larger than 
—dN. 



3 Multicanonical Sampling 

In the pedagogical review H I emphasized that the inverse spectral density is the appropriate 
weight factor to obtain a flat energy density 

w(E) = n-\E) = e -m)E+ a (E)_ (7) 

Here /3(E), a(E) is the multicanonical parameterization 0, Its rationale is to relate 
to the temperature. It should be noted that MC calculations are insensitive to an over-all 
independent factor, i.e. against replacing w(E) by cw(E). In the following I will exploit this 
property from time to time, and not trace back the corresponding multiplicative or additive 
constants. If necessary they may be obtained by introducing a convenient normalization. 
The spectral density may be written as 

n(E) = e s{E \ (8) 
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where S(E) is the microcanonical entropy P5| . The thermodynamical relation for the inverse 
temperature ((3 = T -1 , where my Boltzmann constant convention is k — 1) is 

<-§■ 

For models with discrete energy values this may be translated into 

= SiE + c) -S iE ) 
e 

where e is the smallest possible energy increment such that n(E + e) and n(E) are both 
non-zero. I.e. typically we have e = 4 for the model of section 2 (special care is needed for 
the IF close to its groundstate). Note that equation (10) is in part convention. Other valid 
options would be (3(E) = [S(E) - S(E - e)]/e or (3(E) = [S(E + e) - S(E - e)]/(2e). For 
consistency with 0, [5| I stay with (10). 

Once (3(E) is given, a(E) may be determined recursively. The equality of e~ s ( E > and 

e - P{E )E+ a {E) implies 

S(E) - S(E -e)= (3(E)E - (3(E - e)(E - e) - a(E) + a(E - e). 

Using (10) to eliminate the term e(3(E — e), we find for a(E) the recursion relation 

a(E -e) = a(E) + [f3(E - e) - (3(E)]E, a(E max ) = 0. (11) 

Here a(i? max ) = is a choice of the over-all multiplicative constant, needed to start off the 
recursion. 

To perform a multicanonical simulation, we do not need to know the exact weight factor 
(7). Instead, a working estimate w(E) of w(E) is sufficient, such that the sampled energy 
histogram H(E) is approximately flat in the desired energy range (1). In the subsequent 
discussion I use the notation n(E), S(E), (3(E), a(E) for estimators of the corresponding 
quantities n(E), S(E), (3(E) and a(E). 

The technical feasibility of multicanonical sampling depends on the existence of efficient 
methods to obtain an acceptable estimate w(E). Computational resources and concurrent 



numerical options determine how much computer time one will consider acceptable for the 
calculation of w(E). Typically it should be less or at most of the order of the CPU time 
spent on the subsequent multicanonical sampling (with then fixed parameters). It seems 
that different workers in the field have tried various approaches. I am only familiar with two 
of them. 

(a) Approaches which work in one or two steps fl], ||], 12|. Employing FSS a reasonable 



good approximation W^(E) is obtained by extrapolation from previously simulated, 
smaller lattices. With w^(E) a first multicanonical simulation is carried out. Its 
results give an improved estimate w^ 2 \E) with which additional simulations may be 
done. This approach works well for first order phase transitions, but failed badly for 
some disordered systems. 

b) Recursive calculations ^(E) — > have been employed. They are subject of the 

following subsection. 

3.1 Recursive multicanonical calculations 

Let H n (E) be the unnormalized histogram obtained from a (short) multicanonical simulation 
with w n (E). At energy values for which H n (E) is reliable , the new estimate is 

v ' H n (E) K ' 

Clearly (12) fails for energy values for which H n (E) = 0, and also values like H n (E) = 1 
or 2 are of course statistically unreliable. Worse, even large values like H n (E) = 10 6 may 
still not give reliable estimates. Namely, situations can be encountered where the integrated 
autocorrelation time is of the same order of magnitude or even larger. Before I come to 
a more thorough discussion of this problem, I would like to discuss two approaches in the 
literature. 

To be definite, let us assume that the starting point for the recursion is 

*P{E) = 1. (13) 



In general this is a reasonable choice, which will allow us to recover the normalization (4) 
when desired. For some practical applications other choices, like a canonical simulation at a 
certain temperature, may be more convenient. 

In the paper by Celik and myself [g] equation (12) was stated in the multicanonical 
notation (7). It reads then (note e = 4 in Ref.|||) 

T + \E) = T(E) + e" 1 \n[H n (E + e)/H n (E)}. (14a) 

The function a n+1 (.E') is then determined by equation (11). In addition to (14a) specific 
rules were given about how to exclude unreliable histogram entries. Namely, 

j ^(E) for E > -^median! f-, Ah \ 
{ P (Aut-off J tor & < ^cut-off 

Here -^median i s t ne median of the n th energy distribution, and E™ ut _ oS < -Emodian i s an energy 
cut-off, such that in simulation n the temperature is kept constant for E < E™ ut _ oS . Further, 
note that the starting condition (13) becomes (3°{E) = 0, a°(E) = 0. 
Lee [T^ states his recursion in two parts: 



and 



S n+ \E) = S n (E)+\nH n (E) for H n {E) > 1, (15a) 



S n+ \E) = S n (E) for H n (E) = 0, (156) 



The first part is obviously equation (12), as follows from w n {E) = exp[—S(E)]. The 
identity [14 1 of (15a) and (14a) follows from (10). Obviously (15a) is a convenient interme- 



diate step to derive (14a). The second part (15b) is a specific prescription about how to 
handle H n (E) = 0. The other unreliable H n (E) are included into the recursion (12). Let us 
note the following: 

(a) Besides from minor notational differences, it is uniquely determined how to handle the 
reliable part of the data. One should note that the equivalent equations (12), (14a) and 
(15a) are all non-local in the sense that ultimately histogram entries over the entire 



sampled range will determine the transitions amplitudes from one energy to the next. 
It may be a little surprising that equation (14a) looks less local than equation (12) or 
(15a). This is entirely irrelevant, because the weight factors are only auxiliary quanti- 
ties to determine (for instance by detailed balance) the decisive transition probabilities 
w[E — > E']. The transition probabilities relate different energies. W = (w[E — > E'\) 
forms a (sparse) matrix, and its eigenvector with eigenvalue one is finally supposed to 
become the spectral density, i.e. determines the weight factors. This diagonalization 
(implicitly carried out by the MC simulation) is a non-local process. This non-locality 
may induce certain instabilities. For instance, if inaccurate weight factors somewhere 
in the energy range create a region of attraction, all CPU time in the next run may be 
wasted on iterating on an irrelevant energy region. 

(b) As the recursion (12) and (15a) stand, the statistical accuracy of estimate n + 1 is 
entirely determined by MC simulation n. With increasing n the covered energy range 
gets larger and larger. One needs longer and longer simulations just to regain the 
previously reached statistical accuracy (on the appropriate energy subrange). It is 
possible, but tedious, to combine the statistics of simulations n, n — 1, 1, 0. The 
gain is not quite as dramatic as one may superficially expect, because simulation n+ 1 
does still explore the entire energy range, just for the purpose to explore an additional 
increment of the desired energy range. 

(c) Note that the median rule of (14b) freezes estimates on some part of the already covered 

energy range. One should improve on it by using subsequent statistics when available. 
In [0 it was suggested to combine the median rule with upper bounds on the energy, 
such that the energy range gets reasonably restricted. However, it is then difficult to 
ensure ergodicity. 

(d) A central difficulty of the recursion is the handling of energy regions for which reliable 

statistical information is not yet available. I elaborate on this now. 



S 



Lee's proposal (15b) looks attractive because of its simplicity. It works for the very 
small systems considered in his paper, but for many realistic situations it will lead to an 
unacceptable slowing down. The reason is that (15b) is equivalent to simulating with a 
constant weight factor (7). Now, at low temperatures one typically encounters 

n(E-e)/n(E) ~V~\ (16) 

Therefore, for a not yet covered energy range E < E one will need of order V attempts just 
to achieve once the transition E — > E — e. 

The rule f3 n+l = (3 n+l (E cnt _ oS ) for E < -E cu t- ff from (14b) achieves a far better perfor- 
mance for this situation. Assume that /3(E) is monotonically increasing towards lower ener- 
gies (exceptions are first order phase transitions). A canonical simulation with /? n+1 (.Edit-off) 
will have its maximum energy density at E — E cut _ of[ , because its first derivative with respect 
to the energy is zero there. The width of its energy distribution is of order y/V. Conse- 
quently, there will be no weight factor problem associated with proceeding towards lower 
energies. In practice one has to use estimators (3 (E). One would like to chose -E cu t- ff 
as low as possible, but one encounters noise problem when the cut-off energy is shifted too 
far towards the edge of the reliably covered energy range. With some experience a good 
"pick" for -Ecut-off can be achieved by just inspecting the function (3 (E). Alternatively, 

12 + 1 71+1 

one may use a fit /5 max from several energy values instead of j3 (E , cut _ fr), or even fit the 

72+1 

continuation of the entire function (3 (E) for E < E cut _ oS (with the penalty of spurious 
instabilities). In any case, in energy regions where (16) holds, one expects a performance 
increase by at least a volume factor over using (15b). On the other hand, it is precisely 
this part of the recursion (14) which required annoying attention by hand. This experience 
can, of course, note rule out the possible existence of some more perfect fitting procedure, 

72 + 1 

to estimate (3 (E) towards lower energies. 

How the recursion (14) slows down with volume depends thus on the details of its imple- 
mentation. Typically, one has to cover a macroscopic energy range, i.e. -E max — -Emm ~ V. 
The optimal slowing down of a single multicanonical simulation on this range is ~ V 2 , cor- 
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responding to a random walk in the energy |3J. Of order V ' 5 simulations are needed to 
iterate from an initial canonical distribution up to covering the entire energy range multi- 
canonically. This leads to an optimal slowing down ~ V 2 ' 5 for the recursion. That this is not 
an overestimate follows from the fact that the slowing down of a multicanonical simulation 
on half the energy range still scales with V 2 , and it still takes of order V ' 5 simulations to 
iterate from half the range to the full range. 



4 Accumulative Recursion 

I now introduce a recursion which calculates (3 (E) on the basis of the statistics accumu- 
lated in all previous runs n, n — 1, 1. For this purpose let us first re-write (14a) as 

T + \E) = e 1 \n[H n (E + e)/H${E)\, (17) 

where 

Hp{E) = H n (E)e- iT{E)e . (18) 

Equation (17) still holds when H n (E) and HJ^(E) are replaced by non-zero linear combina- 
tions H n (E) and H$(E): 

n 

H n (E) = W m (E)H m (E), (19a) 

n 

Hp{E) = J2 W m {E)H™{E). (196) 

m=0 

The accumulated statistics can be presented by suitable choice of the weight factors W m (E). 
The optimal choice is not clear, as it may depend non-trivially on the dynamics. In practice 

W m (FV _ min[H m (E + e), H m (E)\ 

K } max[H™(E + e),H™(E)} 1 ; 

has worked well. It relies on the conservative assumption that each contribution to the 
estimate 

P n+ \E) = e- 1 \n[H n {E + e)/H${E)] (21) 
10 



will be as good as its weakest part. This equation is supplemented by 

T + \E) = T + \E + e) (22) 

for the case that either H n (E + e) or Hp(E) has insufficient statistics. To provide some 
feeling for the estimator (21) let me discuss two special cases. 

(a) When the desired, flat distribution is already reached, the weight factors (20) equal 

1 up to statistical fluctuations. Let us ignore fluctuations for the moment. Then 
H n ~ l (E + e) = H n ~ l (E) holds before the n th run, which uses (3 n (E) as defined by 
equation (21). In the n th recursion H n (E+e) = H n (E) is obtained by assumption. This 
leads to H n (E + e) = H 71 ' 1 + H n (E + e) and Hp(E) = H^~\E) + H n (E) exp(-/Te). 
Equations (19), (21) yield p (E) = ~j3 (E), i.e. the /3(E) function is a fixed point 
when the sampled distribution is flat. 

(b) Consider the first recursion, carried out with (3°(E) = 0. The sampling results will be 
H°(E + e)/H°(E) = n(E + e)/n(E), again up to statistical fluctuations. Recursion 
(21) yields f3 1 (E) = e _1 \n[n(E + e)/n(E)], which is already the final multicanonical 
answer due to the fact that we have neglected statistical fluctuations. Quite generally 
it can be shown that the desired multicanonical function (3(E) is an attractive fixed 
point of the recursion. 

In practice there may be severe statistical fluctuations due to only few, correlated en- 
tries in H n (E + e), H n (E) or both. If the number of entries in both arrays is small, but 
approximately equal (W n (E) w 1), equations (19) guarantee that increase from H n ~ l — > H n 
is in proportion the the generated statistics (assuming similar autocorrelation time in runs 
n — 1, n — 2, ...). If the number of entries is only small in either H n (E + e) or H n (E), the 
weight factor (20) correct for the asymmetry. The larger statistics is reduced to the smaller 
one, and the smaller even more suppressed. As the ratio H n (E + e)/H(E) determines the 
estimate (3 n (E) ) it is clear that a large statistical fluctuations in either the numerator or the 
denominator is sufficient to destroy the entire estimate. The weight factor prevents this. 
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The obvious advantage of equation (21) over recursion versions of section 3 is that the 
accumulative statistics of all runs is used to reduce statistical fluctuations. In [|I7[] we have 
not supplemented the present recursion by a median restrictions of the type (14b), although 
this might lead to further improvements. Without such restrictions, typically the recursion 
leads quickly to rather high (3 values, and works its way back from the corresponding low 
energy values through the entire energy range. Occasionally this has led to "hang-up" situ- 
ations, for which a simple "retreat" strategy has turned out to be sufficient. For the case of 
generalized Ising model, the appendix gives and explains an actual program listing, which 
was used for the numerical illustrations of the next section. A generalization of my recursion 
to non-flat distributions, like for instance those proposed in |7j would be straightforward. 



5 Numerical Tests 

I confine myself to reporting results for the 3d IF and the 3d EAI spin glass. Similar tests 
have been performed for the 2d IF and are in progress for the 2d EAI spin glass as well 



as for fully frustrated Ising models p9[ . To keep the relation to the program listing in the 
appendix close, I shall use 

I A = ^(-E + dN), (with iV = L d ) (23) 

instead of the energy, defined by (2). The rationale of I a is its range: 

I a = 0, 1, 2, dN/2 (24) 

in typical increments of 1. For comparison, we had — dN < E < dN in typical increments 
of 4. Consequently, for the purposes of programing I a is far more convenient. Functions of 
E are now interpreted as functions of Ia in the obvious way, i.e. (3(E) = (3[E(Ia)] — > (3(1 a), 
and so on. 
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5.1 Three dimensional Ising ferromagnet 

The first few terms of the low temperature expansion on a finite (but sufficiently large) lattice 
collected in table 1. The present computer program is unsuitable to cope with ti(Ia) = 
for I A = (3N/2) - 1, (3N/2) - 2 and (3N/2) - 4. I just bypass 28 the problem by restricting 
the updating to the range 1a < N max = (3N/2) — 5. Proposals with 1a > N max are simply 
rejected. 

We want to calculate multicanonical parameters for the temperature range infinity down 
to zero. Simulations with (3 = are peaked around I a = N min = 3N/A. We therefore fix the 
function (3 to (3(1 a) = for I a < N min , and never change it there. For I a > N min we perform 
the multicanonical recursion of section 4. The covered range of lattices was 4 < L < 16. In 
a first set of runs the recursion was applied until the system tunneled at least 60 times. The 
(expected) experience from these runs is that the recursion remained stable after the first 
tunneling. The tunneling time tau is then measured after the first tunneling has occurred, 
while continuing to update the parameters. Table 2 collects the measured tunneling times 
r, and states on how many tunneling events n T the estimates actually rely. 

By r I denote the time (as always in updates) it takes until the first tunneling has taken 
place. This is essentially the time our recursion needs to provide a reliable estimate of the 
multicanonical parameters, and it will therefore be called recursion time in the following. 
Two estimates, Tq and Tq, are given in table 2. They differ by the number of sweeps 
performed before the multicanonical parameters are updated (i.e. the subroutine UPMUCA 
of the appendix is called). A sweep is defined by updating N = L d spins. For Tq UPMUCA 
was called every 120 sweeps, whereas for r b it was called every iV sweeps. Respectively, this 
amounts to letting the numbers of intermediate updates grow in proportion to the volume 
and to the volume squared. Within the (still large) statistical errors there is no difference 
noticeable. 

The values n T a and n T b are the numbers of (3(1 a) = re-starts on which the respective 
estimates rely. As the average CPU time needed per recursion is substantially higher than 
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the average tunneling time r, I have limited the tq analysis to L < 12. The given error 
bars are somewhat unreliable as the obtained distributions have long tails towards large To 
values. A large statistics is needed to get into the region where the central limit theorem 
provides a good approximation. My typical number of n T0 = 126 events is a bit at the low 
edge. Figure 1 employs a log-scale for r to show the histograms for Tq. The distributions 
for the tunneling times r themselves, are more reasonably, Poisson like, behaved. 

Figure 2 shows the increase of r and Tq with volume on a log-log scale. The straight 
lines correspond to the fits r = cV 5 and Tq = coV s °. The results for the fit parameters are 

ln(c) = -0.53 ±0.16, 6 = 2.249 ± 0.021, (Q = 0.18) (25) 

and 

ln(c ) = -1.24 ± 0.17, 5 = 2.931 ± 0.023, (Q = 0.70), (26) 

where Q is the goodness of fit @]. It should be remembered that the lower bounds are 6 = 2 
[H] and 5o = 2.5 (see section 3). 

To demonstrate that after a few tunneling events the multicanonical parameters are 
indeed already useful, I have also measured a tunneling time T\, obtained by fixing the 
multicanonical parameters after the first four tunneling events. Table 2 contain also the 
corresponding estimates T\. Within the statistical errors, there is no difference with the 
estimates r. 

5.2 Three dimensional Edwards— Anderson Ising spin glass 



First I present some results from the extensive investigation [|17|], which are not contained 
in this reference. For fixed lattice size L tunneling times r are found to vary greatly for 
different realizations. For each lattice size figure 3 connects the tunneling times, sorted 
in decreasing order. For L = 4 — 8 there are 512 different realizations per lattice. For 
L = 12 there are only seven realizations, depicted at 64(i — 1), (i = 1, 7). The lines are 
drawn to guide the eyes. Figure 4 depicts histograms for the L = 4 — 8 tunneling times. 
In both figures a logarithmic scale is used for r. The worst realizations have dramatically 

14 



larger tunneling times than typical ones, defined by the median value T0.5. This leads to 
large differences between the mean value r, which determines the needed computer time, 
and the median value T0.5. These values are collected in table 3. With increasing lattice size 
the discrepancy between mean and median increases dramatically (the L = 12 data have to 
be considered unreliable for this purpose). This lack of self-averaging of the spin glass with 
respect to the multicanonical tunneling time comes somewhat surprising, and remains to be 
better understood. Also collected in the table are the smallest r .o and largest r L0 tunneling 
time, found on the investigated realizations. 

For typical spin glass realizations, i.e. the realizations corresponding to the median t 05 
tunneling times of table 3, I have performed the same analysis as for the 3d IF in the previous 
subsection. The results are collected in table 4. An interesting and unexpected result is that I 
find T\ systematically smaller than r, i.e. further applications of the recursion relation make 
the tunneling worse. My tentative interpretation is that the flat distribution is not optimal. 
Due to statistical fluctuations, one can then imagine that immediately after one of the 
first few tunneling events the generated multicanonical parameters are positively correlated 
towards a more optimal choice. A more detailed future analysis may be desireable. 

As before, the recursion times Tq and t\ are practically identical. However, a second 
unexpected result is that now the recursion times take the same order of magnitude as the 
tunneling times, whereas for the 3d IF the recursion times were considerably large than the 
tunneling times. 

It has to be remarked that the L — 12 results are not in line with the subsequent 
estimates from lattices of size L = 4, 6 and 8. (a) The f\ estimate is considerably higher 
than expected. The reason is likely that the typical realization picked is not typical. A 
reliable estimate of To. 5 is practically impossible due to the small number of only seven 
L = 12 realizations investigated in |T7[ . (b) The Tq value, given in brackets, is much smaller 
than expected. However, the number is given in brackets as it is not an estimate of said 
quantity. Altogether twelve attempts were made, to determine multicanonical parameters 
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by means of the recursion. Of those, two did not lead to a single tunneling event within the 
maximally allowed CPU time corresponding to approximately 11.2E08 updates. These two 
attempts are (cannot be) not included in the given average values. This behavior illustrates 
that one should perform several independent starts, when applying the recursion to difficult 
situations. 

Using only the estimates from L = 4, 6 and 8, the subsequent my results are obtained 
from straight line fits to the equations r = cV 5 , t\ = c\V Sl and Tq = cqV 5 °: 

ln(c) = -3.04 ±0.29, 5 = 3.24 ± 0.06, (Q = 0.40), (27a) 

ln(ci) = -3.61 ±0.24, 5 = 3.28 ± 0.05, (Q = 0.39), (276) 

and 

ln(co) = -2.23 ± 0.24, 5 = 3.09 ± 0.04, (Q = 0.78), (28) 

Here, as well as in the previous section, the routine GFIT from gives results perfectly 
compatible with the linear fit results. A figure corresponding to (27) and (28) looks similar 
to figure 2, but is not very instructive as all three fits lines are almost on top of one another. 



The exponent 5 is smaller than the one reported in ||17|| . The reason is that it is differently 
defined. In |T7[] the tunneling time was averaged over all realization, whereas here I have 
picked single, typical realizations. There is evidence that for the worst realizations the tun- 
neling time slows down exponentially with L. This spoils the power law fit for the average 
over all realizations. 



6 Summary and Conclusions 

For the 3d Ising ferromagnet it is clear that the FSS methods employed in [|l], [| provide 
reliable estimates of the multicanonical parameters more efficiently than the recursion of 
this paper. On the other hand, the FSS approach breaks down || for the important class 
of disordered systems. Then recursions like the one of this paper become crucial to enable 
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the method, and the Ising ferromagnet is still a suitable testing ground to set quantitative 
performance scales. These are now given, for the first time, by tables 2 and 4. Table 4 
corresponds to the important case of a typical Edwards-Anderson Ising spin glass. Future 
investigations will have to cope with these standards. It is my hope that they will bring 
improvements in the constant factor, and possibly towards a V 2 power law behavior, which 
is optimal for any kind of local random walk behavior. 



7 Appendix 



In this appendix I describe the actually used computer implementation for the accumulative 
recursion of the multicanonical parameters. The relevant Fortran subroutine is listed next. 
It is not claimed that this subroutine is in any sense optimal. It just worked sufficiently well 
for the described examples. 

SUBROUTINE UPMUCA(IRPT) 
C Update of multicanonical parameters. 
C HAMUA(*,1): over-all sum (record keeper only). 
C HAMUA(*,2): LRTRT adjusted over-all sum (record keeper only). 
C HAMUA(*,3): 1. weighted sum. 
C HAMUA(*,4): 2. weighted sum. 

IMPLICIT REAL* 8 (A-H.O-Z) 

IMPLICIT LOGICAL (L) 

PARAMETER (ND=3 , NL=08 , NS=NL**ND , NRPT=100 , NSW=NS) 

PARAMETER (NNH= (ND*NS) /2 , NAMIN=NNH/2 , FRTRT=3 . DO , EPS=1 . D-8) 

PARAMETER (HMIN=1 . 0D00*FL0AT (NS) *FL0AT (NSW) ) 

COMMON /MEAH/ HA (0 : NNH) , IAMIN , IAMAX , ITMIN , ITMAX 

COMMON /MUCA/ B (0 : NNH) , A(0 : NNH) ,HAMU(0 : NNH ,4) , LRTRT (NRPT) 

C 

DO IA=ITMIN, ITMAX 
HAMU(IA,1)=HAMU(IA,1)+HA(IA) 
HAMU ( I A , 2 ) =HAMU ( I A , 2 ) +HA ( I A) 
END DO 

C 

C Retreat strategy (below) implies: range up to IAMAX . GE . ITMAX . 
IAMAM1=IAMAX-1 
DO IA=NAMIN , IAMAM1 
IAP1=IA+1 

HAMIN=MIN (HA (IA) ,HA(IAP1)) 
H AMAX=MAX ( H A ( I A ) ,HA(IAP1)) 

IF(HAMIN.GT.0.5D00) THEN 
W1=HAMIN/HAMAX 

HAMU(IA,3)=HAMU(IA,3)+W1*HA(IA) 
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HAMU(IA,4)=HAMU(IA,4)+W1*HA(IAP1)*EXP(-4.0D00*B(IAP1)) 

END IF 

C BETA update (after retreat HAMIN.LE.0.5 possible): 
HAMUMIN=MIN (HAMU ( I A , 3) , HAMU ( I A , 4) ) 

IF (HAMUMIN.GT.EPS) THEN 
B(IAP1)=-0.25D00*L0G(HAMU(IA,4)/HAMU(IA,3)) 

ELSE 

B(IAP1)=B(IA) 

END IF 

END DO 

C 

C Retreat strategy for hung-up situations: 

LRTRT(IRPT)=. FALSE. 
C Besides retreat, update of MUCA A-array is performed 
C (range up to IAMAX.GE. ITMAX is needed for this reason). 

DO IA=NAMIN , IAMAM1 

IAP1=IA+1 

IF(HAMU(IAP1,2) . GT . HMIN . AND . HAMU ( I AP 1 , 2 ) . GT . FRTRT*HAMU( IA , 2) ) THEN 
C The program may need modifications, if there are 
C energy values without states in the .LE.IAMAX range. 

IF(HAMU(IA,2) .EQ.O) PRINT* , 'UPMUCA Warning: IA = ',IA 

IF( .NOT . LRTRT(IRPT) ) PRINT*, 
& 'RETREAT! IRPT , IA , HAMUs : ' , IRPT, IA ,HAMU(IAP1 , 2) ,HAMU(IA,2) 

LRTRT ( IRPT) = . TRUE . 

END IF 

IF (LRTRT (IRPT)) THEN 

HAMU ( I AP 1 , 2 ) =HAMU ( I AP 1 , 2 ) /FRTRT 

B(IAP1)=0.0D00 

HAMU ( I AP 1 , 3 ) =HAMU ( I AP 1 , 3 ) /FRTRT 
HAMU ( I AP 1 , 4) =HAMU ( I AP 1 , 4) /FRTRT 

END IF 

A(IAP1)=A(IA)-4.0D00*(B(IAP1)-B(IA))*FL0AT(IA) 
END DO 

C 

IAMAP1=IAMAX+1 
DO IA=IAMAP1,NNH 
B(IA)=B(IA-1) 
A(IA)=A(IA-1) 
END DO 

C 

RETURN 
END 



Relevant parameters (to be set) are the dimension ND, and lattice size NL. The presented 
choice is an 8 3 lattice. NS encodes the lattice size and NNH is needed to dimension a number 
of arrays. 

The argument IRPT keeps track of the number of repeated calls to UPMUCA. In an outside 
DO-loop IRPT runs from 1 to NRPT. Inside our subroutine NRPT is only needed to dimension 
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the LOGICAL array LRTRT, which keeps track of the number of "retreats", to be discussed 
later. A parameter not needed at all in our subroutine is NSW. It denotes the number of 
update sweeps performed in between the calls to UPMUCA. In the presented code it is set 
equal to the lattice size, corresponding the recursion time Tq of section 5. 

NAMIN sets the lower bound on the IA range (I a of section 5) to which the recursion is 
applied: B(IA)=0 for IA < MAM IN implements 13(1 a) = for I a < N min . The other parameters 
will be discussed later on. 

Most arguments are passed through COMMON blocks. On entry the array HA contains 
the newly assembled statistics, i.e. the histogram of the number of times a certain IA 
value (corresponding to an energy via (23)) has been visited during the last NSW sweeps. 
(The information is collected after each single spin update. The array HA has to be set to 
HA = after each call to UPMUCA.) Further arguments passed by the COMMON block MEAH 
(measurements) are: IAMIN, the smallest IA value encountered so far (not used in UPMUCA); 
IAMAX, the largest IA value encountered so far; ITMIN, the smallest IA value encountered 
during the last NSW sweeps, and ITMAX, the largest IA value encountered during the last NSW 
sweeps. 

The meaning of the array(s) HAMU is explained by the comments at the beginning of the 

subroutine. Central for the code are the lines 
W1=HAMIN/HAMAX 

HAMU(IA,3)=HAMU(IA,3)+W1*HA(IA) 

HAMU(IA,4)=HAMU(IA,4)+W1*HA(IAP1)*EXP(-4.0D00*B(IAP1)) 

which implement our equations (19) and (20) recursively. Next, the arrays A and B 
correspond to the multicanonical functions (5 and a the lines 

B(IAP1)=-0.25D00*L0G(HAMU(IA,4)/HAMU(IA,3)) 

and 

A(IAP1)=A(IA)-4.0D00*(B(IAP1)-B(IA))*FL0AT(IA) 

implement equations (21) and (11). Of course, A(NAMIN)= 0. The parameter EPS pre- 
vents that the /^-recursion takes place without sufficient statistics, and otherwise equation 
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(22) is chosen. 

Some complications arise, mainly because a "retreat" strategy has been implemented to 
get out of certain "hung-up" situations. To discuss them is beyond the scope of this paper, 
as the relevant (spin glass) configurations require more detailed investigations first. In short, 
an extreme difference between HAMU(IA+1,2) and HAMU(IA,2) can turn out to be artificial, 
such that its statistics is better not trusted. "Extreme" is defined by the parameter FRTRT, 
put to 3 in the presented code. When the thus defined limit is exceeded the assembled statis- 
tics is reduced in weight by the factor 1/FRTRT and (3(1 a) is put in the corresponding energy 
region to 0(1 a) = for the next recursion. The program may thus escape certain traps 
successfully. However, I like to remark that one has to chose FRTRT to be very large (around 
200), if one likes to calculate multicanonical parameters for an 24 3 IF in the range described 
in section 5.1. The reason is the peculiar IF density of states anomaly from I A = (3N/2) — 7 
to 1a = (3N/2) — 6 (see table 1). The choice (3 = re-creates an FRTRT factor of order N. 
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E 


Ia 










-3N 


3iV/2 


2 


-3iV + 4 


(3JV/2) - 1 





-3iV + 8 


(3JV/2) - 2 





-3N + 12 


(37V/2) - 3 


2N 


-3N + 16 


(3JV/2) - 4 





-3N + 20 


(3iV/2) - 5 


6iV 


-3iV + 24 


(3JV/2) - 6 


2N 2 - UN 


-3N + 28 


(3iV/2) - 7 


30iV 


-3N + 32 


(3JV/2) - 8 


6iV 2 - 66iV 


-3iV + 36 


(3JV/2) - 9 


(27V 3 -42iV 2 + 1252iV)/6 



Table 1: Finite lattice low temperature expansion for the 3d IF (N = L 3 ). 
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T a 
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4 


548 


719 (19) E01 


126 


661 (41) E02 


126 


557 (46) E02 


111 


731 (58) E01 


6 


354 


126 (05) E03 


252 


195 (20) E04 


126 


219 (23) E04 


145 


129 (10) E03 


8 


559 


881 (23) E03 


126 


311 (55) E05 


126 


253 (50) E05 


125 


839 (69) E03 


12 


322 


118 (06) E05 


140 


95 (32) E07 


164 


90 (15) E07 


141 


127 (11) E04 


16 


577 


760 (30) E05 


180 




2 


14 (big) E09 


180 


746 (54) E05 



Table 2: Tunneling and recursion times for the 3d IF. 
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T 0.0 


T 0.5 


Tl.O 












4 


398 (15) E02 


144 E02 


304 E02 


411 E03 


6 


336 (30) E04 


436 E03 


131 E04 


670 E05 


8 


171 (46) E06 


505 E04 


282 E05 


213 E08 


12 


139 (77) E08 


408 E06 


481 E07 


544 E08 



Table 3: Mean r and some q-tiles r q for the 3d EAI tunneling time. 
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4 


185 


332 (25) E02 


126 


523 (19) E02 


126 


409 (22) E02 


270 


228 (14) E02 


6 


256 


181 (12) E04 


126 


150 (10) E04 


126 


172 (11) E04 


357 


131 (08) E04 


8 


134 


272 (26) E05 


252 


245 (13) E05 


252 


253 (16) E05 


207 


203 (16) E05 


12 










10 


[19 (4) E06] 


13 


92 (31) E08 



Table 4: Tunneling and recursion times for typical 3d EAI spin glass realizations. 
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Figure 1: Histograms for the 3d IF recursion time Tq on lattices of size L 3 . 
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This figure "figl-2.png" is available in "png" format from: 



http://arXiv.org/ps/hep-lat/9503019v2 
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ure 3: Sorted 3d EAI tunneling times for various lattices of size L 3 . 
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Tunneling Time Estimates: Histograms over Realizations. 
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Figure 4: Histograms over realizations for 3d EAI tunneling time estimates 
T on lattices of size L 3 . 
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